# process all images itarable
options(warn=-1)
library("raster")
library("sf")
library("tidyverse")
library("rasterVis")


config.ndwi_treshold = 0 # allowed: -0.05, 0, 0.1 by article 
config.mask <- st_read("data/mask/tile.geojson")

# list directories by years
years = list.dirs("data/sattelite", recursive = F, full.names = F)

df.year <- c()
df.month <- c()
df.day <- c()
df.ndci_mean <- c()
df.ndci_min <- c()
df.ndci_max <- c()
df.chl_mean <- c()
df.chl_min <- c()
df.chl_max <- c()

# iterate years
for (year in years) {
  print(paste0("======== Processing year: ", year))
  mds <- list.dirs(paste0("data/sattelite/", year), recursive = F, full.names = F)
  # iterate month-day folders
  for (md in mds) {
    mdsp <- str_split(md, '-', simplify = T)
    month <- mdsp[1]
    day <- mdsp[2]
    print(paste0("= Processing day: ", month, "-", day))
    
    # iterate files in folder
    files = list.files(paste0("data/sattelite/", year, "/", md))
    fb3 <- NULL
    fb4 <- NULL
    fb5 <- NULL
    fb6 <- NULL
    for (f in files) {
      if (grepl("B03_", f)) {
        fb3 <- paste0("data/sattelite/", year, "/", md, "/", f)
      }
      if (grepl("B04_", f)) {
        fb4 <- paste0("data/sattelite/", year, "/", md, "/", f)
      }
      if (grepl("B05_", f)) {
        fb5 <- paste0("data/sattelite/", year, "/", md, "/", f)
      }
      if (grepl("B08_", f)) {
        fb8<- paste0("data/sattelite/", year, "/", md, "/", f)
      }
    }
    
    if (is_null(fb3) || is_null(fb4) || is_null(fb5) || is_null(fb8)) {
      print(paste0("Raster files are missing: ", year, '-', md))
      next
    }
    
    b3 <- raster(fb3)
    b4 <- raster(fb4)
    b5 <- raster(fb5)
    b8 <- raster(fb8)
    
    b3 <- raster::mask(b3, config.mask)
    b4 <- raster::mask(b4, config.mask)
    b5 <- raster::mask(b5, config.mask)
    b8 <- raster::mask(b8, config.mask)
    
    ndwi <- (b3 - b8)/(b3 + b8)
    # make binary reclassification to 2 classes: from -1 to 0.3 => 0, from 0.3 to 1 => 1
    cndwi <- reclassify(ndwi, c(-1, config.ndwi_treshold, 0, config.ndwi_treshold, 1, 1))
    
    mask_ndwi <- rasterToPolygons(cndwi, fun = function(x){x == 1}, digits=6, dissolve = T)
    shapefile(mask_ndwi, paste0("output/rasters/ndwi/", year, "/", month , "-", day , ".shp"), overwrite=T)

    b3 <- raster::mask(b3, mask_ndwi)
    b4 <- raster::mask(b4, mask_ndwi)
    b5 <- raster::mask(b5, mask_ndwi)
    b8 <- raster::mask(b8, mask_ndwi)
    
    ### NDCI CALC
    # calc chlorophyll index (normalized)
    ndci <- (b5-b4)/(b5 + b4)
    # make custom extent
    e <- extent(44.2, 44.7, 45.45, 45.60)
    ndci <- crop(ndci, e)
    writeRaster(ndci, paste0("output/rasters/ndci/", year, "/", month , "-", day , ".tiff"), overwrite=T)
    
    # plot output NDCI
    png(paste0("output/maps/ndci/",year,"/",month,"-",day,".png"), res=180, width = 1200, height=850)
    pal <- terrain.colors(4)
    plot(ndci, breaks = c(0, 0.2, 0.4, 1), col=pal)
    dev.off()
    # process values
    vals.i <- as.vector(na.omit(ndci@data@values))
    stats.ndci = list(
      mean = mean(vals.i),
      sd = sd(vals.i),
      stderr = sd(vals.i)/sqrt(length(vals.i)),
      min = min(vals.i),
      max = max(vals.i)
    )
    png(paste0("output/hist/ndci/",year,"/",month,"-",day,".png"), res=180, width = 1200, height=850)
    hist(vals.i)
    dev.off()
    
    ### Calc Chlorophyll-a by Xing Wang et al., 2018
    ugl <- 20.877*ndci + 3.273
    writeRaster(ugl, paste0("output/rasters/chla/", year, "/", month , "-", day , ".tiff"), overwrite = T)
    png(paste0("output/maps/chla/",year,"/",month,"-",day,".png"), res=180, width = 1200, height=850)
    pal <- terrain.colors(5)
    ugl@data@values[ugl@data@values < 0] <- 0
    plot(ugl, col=terrain.colors(15), main=paste0("Хлорофилл-а ", year, "-", month, "-", day))
    dev.off()
    vals <- as.vector(na.omit(ugl@data@values))
    vals[vals < 0] <- 0
    stats.chl = list(
      mean = mean(vals),
      sd = sd(vals),
      stderr = sd(vals)/sqrt(length(vals)),
      min = min(vals),
      max = max(vals)
    )
    
    png(paste0("output/hist/chla/",year,"/",month,"-",day,".png"), res=180, width = 1200, height=850)
    hist(vals)
    dev.off()
    
    df.year <- c(df.year, year)
    df.month <- c(df.month, month)
    df.day <- c(df.day, day)
    df.ndci_mean <- c(df.ndci_mean, stats.ndci$mean)
    df.ndci_min <- c(df.ndci_min, stats.ndci$min)
    df.ndci_max <- c(df.ndci_max, stats.ndci$max)
    df.chl_mean <- c(df.chl_mean, stats.chl$mean)
    df.chl_min <- c(df.chl_min, stats.chl$min)
    df.chl_max <- c(df.chl_max, stats.chl$max)
    
  }
}

df <- data.frame(df.year, df.month, df.day, df.ndci_min, df.ndci_mean, df.ndci_max, df.chl_min, df.chl_mean, df.chl_max)
write.csv(df, "output/summary.csv", row.names = F)
